submitted to: Japan J. Indust. Appl. Math. 
running title: Curvature adjusted tangential velocity 



Evolution of plane curves with a curvature adjusted tangential velocity* 



Abstract. We study evolution of a closed embedded plane curve with the normal 
velocity depending on the curvature, the orientation and the position of the curve. 
We propose a new method of tangential redistribution of points by curvature adjusted 
control in the tangential motion of evolving curves. The tangential velocity distributes 
grid points along the curve not only uniform but also lead to a suitable concentration 
and/or dispersion depending on the curvature. Our study is based on solutions to 
the governing system of nonlinear parabolic equations for the position vector, tangent 
angle and curvature of a curve. We furthermore present a semi-implicit numerical 
discretization scheme based on the flowing finite volume method. Several numerical 
examples illustrating capability of the new tangential redistribution method are also 
presented in this paper. 

Key Words: curvature driven flow of a plane curve, curvature adjusted tangential 
velocity, semi-implicit scheme, flowing finite volume method, crystalline curvature flow 
equation. 
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1 Introduction 

The purpose of this paper is to study evolution of a plane curve r*,t e [0,T), driven by 
the normal velocity v which is assumed to be a function of the curvature k, tangential 
angle u and position vector a; G T*, 



Geometric equations of the form (11. ip often arise from various applied problems like e.g. 
the material science, dynamics of phase boundaries in thermomechanics, computational 
geometry, semiconductors industry, image processing and computer vision, etc. For a 
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comprehensive overview of industrial applications of the geometric equation having the 
form of (11. ip we refer to a book by Sethian |26j . 

We consider a simple, embedded and closed plane curve F which is parameterized by 
a smooth function x : M/Z D [0, 1] — )■ such that T = Image(a;) = {x(u)] u G [0, 1]} 
and \dux\ > 0. Hereafter, we denote (9^F = df/d^, and \a\ = ■ a where a ■ 6 is the 
Euclidean inner product between vectors a and h. The unit tangent vector can be defined 
as t = dux/\dux\ = dgX, where s is the arc-length parameter and ds = \dux\du. The unit 
inward normal vector n is defined in such a way that det(t, n) = 1. The signed curvature 
in the direction n is denoted by k. With this orientation we have k = det{dsX,dgX). 
Let u be the angle of t, i.e., t = (cos z/, sin i/)""" and n = (— sin i/, cos i/)"'". The problem 
of evolution of curves is stated as follows: For a given initial curve F° = Image(a3°), 
find a family of planar curves {F*}o<t<T, F* = {x{u,t); u G [0, 1] C M/Z} starting from 
F° = {x{u,0) = x^{u)] u G [0, 1]} and evolving in the normal direction according to the 
velocity given by (11.11) . 

As typical examples of /3 one can consider /3 = w{v)k representing the anisotropic 
mean curvature flow, /3 = \k\"'~^k (power like flow arising in image processing) and /3 = 
w{x, i>)k + F{x, u), where w is a weight function, F is an external force and 7 > is a real 
parameter. The normal velocity v is the normal component of the following evolutionary 
equation for the position vector x: 

dtx = (3n + at, x{-,0) = x^{-). (1.2) 

Here a is the tangential component of the velocity vector. Notice that the motion in the 
tangential direction with a tangential velocity a has no effect of the shape of evolving 
closed curves [7,, Proposition 2.4]. The shape is determined by the value of the normal 
velocity /3 only. Hence the trivial setting a = can be chosen. In [6] Dziuk investigated 
a numerical scheme for (5 = k and a = 0. However, it was documented by many authors 
(see e.g. [TOl HH [151 129] and references therein) that such a choice of a may lead to 
various numerical instabilities caused by either undesirable concentration and/or extreme 
dispersion of numerical grid points. 

In order to construct a stable numerical computational scheme, various choices of 
a nontrivial tangential velocity have been proposed and analyzed by many authors. We 
present a brief review of development of nontrivial tangential velocities. In [121 [13] Kimura 
proposed a uniform redistribution scheme for the case when /3 = by using a special 
choice of a satisfying the uniform redistribution condition (U): \dux\ = L*, where L* is 
total length of a curve F*. The author also proved convergence of a numerically discretized 
solution when assuming uniform distribution of initial grid points. In |10], Hou, Lowen- 
grub and Shelley (see also Hou, Klapper and Si [11]) utilized the condition (U) directly for 
the case (5 = k. More precisely, they derived the form of a tangential velocity with ip = 1 
and = (see our notation (13. 3p described in Section [3]). Such a tangential velocity 
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was also proposed and analyzed independently by Mikula and the first author in 
In [lOl Appendix 2], Hon et al. also pointed out the role of tangential velocity (13. 3p for 
a general class of the curvature shape function <f{k) with w = as a generalization of 
condition (U). However, there was no systematic explanation of the consequence of such 
a choice of the tangential velocity on grid points redistribution. In the present paper, the 
important role of the so-called curvature adjusted tangential velocity fl3.3p is emphasized 
from qualitative and quantitative point of view. 

As we already mentioned, in the paper [15j , the authors derived the expression for the 
tangential velocity f l3.3p with ip = 1 and u = for a rather general class of planar mean 
curvature flows satisfying the geometric equation (11. ip with /3 = k). This result can 
be considered as an improvement of that of [H] in which satisfactory results were obtained 
only for the case when /3 = /3{k) is linear or sublinear function with respect to k. Next, 
in the series of papers [161 US [H], Mikula and the first author proposed a method of 
asymptotically uniform redistribution. In terms of our notation, they derived (13.30 with 
if = 1 and nontrivial relaxation function uj{t) for a general class of normal velocities of 
the form (3 = f3{x, z/, k). Their method was also applied to geodesic curvature flows and 
image segmentation problems. 

Besides these uniform or asymptotically uniform redistribution methods, in the so- 
called crystalline curvature flow, the grid points are distributed densely (sparsely) in those 
part of the curve where the absolute value of the curvature is large (small). Although this 
redistribution is far from being uniform, numerical computation is quite stable. One of the 
reasons for such a behavior is that polygonal curves are restricted to a class of admissible 
facet directions. In order to extract essence of the crystalline curvature flow of polygonal 
curves and generalize it to a wide class of plane curve evolution, the second author showed 
that the tangential velocity a = —dsP/k is implicitly involved in the crystalline curvature 
flow of planar curves (cf. [29] )• Notice that such a tangential velocity is a local one since 
its value at some point of a curve depends on the local properties of the curve near x. 
In this point, it is worthwhile to mention the paper by Deckelnick [5] who used locally 
defined tangential velocity having the form a = —du{\dux\~^) for the special case when 
(3 = k. Recently, Pans et al. successfully apphed curvature adjusted tangential velocity 
for evolution of open curves having important applications in the dislocation dynamics 
(see [221123]). 

The asymptotically uniform redistribution is quite effective and applicable for a wide 
range of applications. However, from the approximation point of view, unless the solution 
curve is a circle, there is no reason to accept uniform redistribution for a general case of 
a curve evolution. Hence the desired redistribution should take into account the shape 
of evolution curves and it should also depend on the modulus of the curvature. We will 
furthermore present a combination of the method of asymptotic uniform redistribution [16] 
and the crystalline tangential velocity discussed in |29j . 
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As far as 3D implementation of tangential redistribution is concerned, in the recent 
paper by Barrett, Garcke and Niirnberg [1], they proposed and studied a new efficient 
numerical scheme for evolution of surfaces driven by the Laplacian of the mean curvature. 
It turns out that a uniform redistribution of tangential velocity vectors is implicitly built 
in their numerical scheme. A similar scheme with an explicit expression for the tangential 
redistribution has been proposed and utilized by Morigi in [20] . 

The organization of the paper is as follows: In Section [2l we recall the system of 
PDEs for geometric equations (11. ip and (II. 2p . In Section [3l we will address the question 
how to control redistribution of grid points by using a nontrivial tangential velocity. We 
focus on a class of tangential velocities dependent locally on the curvature combined 
with known asymptotically uniform redistributions. The aim of Section H] is to further 
motivate our study by constructing a curvature adjusted redistribution yielding the best 
possible approximation of the evolving family of planar curves as far as the minimization 
of the error between the length of a curve and the length and enclosed of its polygonal 
approximation is concerned. In Section [5l a numerical solution to the system of governing 
equations will be constructed by means of the so-called flowing finite volume method. 
In Section |6l we will present various numerical examples of applications of the curvature 
adjusted tangential velocity. Finally, we conclude with some key points in Section [71 

2 Governing equations 

Without loss of generality, we will rewrite the normal velocity P{x, u, k) as follows: 



The tangential velocity a will be defined later in Section [31 Using Frenet formula d^x = 
dgt = kn, the following equation for the position vector follows from (11. ip and (II. 2p : 



where the operators dg and stand for arc-length derivatives, i.e. dsF = g~^duf and 
9gF = g'^^dudgf, respectively. Here g = \dux\ denotes the local length of a curve param- 
eterized by X. 

Remark 2.1 If /3 is a linear or superlinear function with respect to k like e.g. w{v)k, 
w{x,i')k + F{x,v), \k\'^~^k (7 > 1), then the function w appearing in (12. ip has no sin- 
gularity at /c = 0. However, in the sublinear case, a singularity occurs at A; = 0. Then a 
regularization at A; = will be necessary from both theoretical (local existence and unique- 
ness of solutions) as well as practical (construction of a stable numerical approximation 
scheme) point of view. We refer the reader to [15] for further discussion of this issue. For 
example, ii w = \k\"'~^^ 7 G (0,1), one can regularize |/c|~^ as follows: \k\~^ = \k\~^ for 
\k\ > e and Ik]"^ = e^^ for \k\ < e, for a small e > 0. 



(3 = w{x, u, k)k + F{x, z/). 



(2.1) 



dtX = wdl 



X + adsX + Fn 



(2.2) 
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Following [15], one can derive a closed system of PDEs governing the motion of curves 
satisfying the geometric equation (11. ip and (11. 2p . These equations can be derived using 
the facts that k = det{dsX,dgX), t = dgX = (cos i/, sin z/)'^, Frenet's formulae dst = kn, 
dgTi = —kt and the relation k = dgV- The resulting system of governing PDEs together 
with (12. 2p reads as follows: 

dtg = {-kl3 + d,a)g, (2.3) 

dtk = dlP + adsk + (2.4) 

dtP = {dkf3)dlp + (a + d,l3)dsP + V./3.t, (2.5) 

for u e [0, 1] and t e (0, T), with initial conditions g{-, 0) = g°{-), k{-, 0) = A;°(-), z^(-, 0) = 
//"(■), x{-,0) = x^{-), and periodic boundary conditions for g,k,x for u G [0,1] C M/Z 
and the boundary condition z/(l,t) = i^(0,t) + 2n. The partial derivative VxP is defined 
as Va;/3 = {dx-i^(3{x, ■, ■),dx2(3{x, ■, ■))'^ for a; = {xi,X2)'^. The initial functions g^,k^,t^ 
and a?" are required to satisfy the following compatibility conditions: g^ = \duX^\ > 0, 
= = {g'>)-^det{dux'^,dlx'>). 

In our simple and fast numerical approximation scheme proposed in Section [5] we will 
discretize (12.21) and the tangential velocity equation (13. 3p with the constraint (13. 4p only. 
With regard to [15j , in the case where a more accurate numerical scheme and results are 
required, it is recommended to discretize the full system of evolutionary PDEs (12. 2p . (12. 3p . 
(12. 4p and (12. 5p for all remaining geometric quantities g, k, and u. On the other hand, as 
far as the time complexity of computation and simplicity of a numerical approximation 
scheme are concerned, discretization of (12. 2p together with (13. 3p . (13. 4p . yields satisfactory 
numerical results. 

However, as w may depend on k and u, the proof of existence and uniqueness of 
a smooth solution to the system of equations (12.20 . (13. 3p . (13.40 does not seem to be a 
straightforward issue. 

Recall that in the case where a is given by the expression (13. 3 p with ip = 1 and u = 0, 
the short time existence and uniqueness of smooth solutions to the system of PDEs (12. 2p . 
(12. 3p . (12. 4p and (12. 5 p has been shown in [15]. In the forthcoming paper [25j the authors 
proved local existence and uniqueness of a classical solution to the full system of governing 
equations ^^-^M and ([S3D, (E3D- 

3 A curvature adjusted tangential redistribution of grid points 

A key tool for construction of a suitable tangential redistribution functional a is the 
so-called relative local length quantity (cf. pL5j). It is defined by the following ratio r: 

riu,t) = ^^, ue[0,l], te[0,T), 
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where L is the total length of F : 

= [ ds= [ gdu, t e [0,T), 
ir* Jo 

and T > is the maximal time of existence of a solution. 

The role of r can be explained as follows: suppose that grid points {x{ui, t)}fLi are 
distributed on the curve F* for Ui = i/N, i = 1, 2, ■ ■ ■ ,N. Since the arc-length s is given 

by s{u, t) = s(0, t) + g{u, t) du, we have 
Jo 

j-Ui j-Ui 

s{ui,t) — s{ui^i,t) = / g{u,t)du = L^ / r{u,t)du 

for each i. Hence, if r{u,t) = 1 holds for all u at a time t, then all grid points are 
distributed uniformly in the sense that s{ui, t) — s{ui-i, t) = L^/N for each i at the time t. 
Furthermore, if the evolving family of curves fulfills the limit r{u, t) — )■ 1 for all m as t — ?■ T 
then redistribution of grid points become asymptotically uniform as t — )■ T (cf. p!6| [T7]). 

It has been documented by many practical examples (see e.g. [HI [17] ) that (asymp- 
totically) uniform redistribution can significantly stabilize and speed up numerical com- 
putation. We recall that the tangential motion has no influence on the shape of evolved 
closed planar curves. A natural question arises: 

Question 3.1 How to design the relative local length ratio r and, subsequently, a such 
that redistribution takes into account the shape of the limiting curve. In other words, 
how to densely (sparsely) redistribute grid points on those parts of a curve where the 
modulus of the curvature is large (small). 

The modulus of curvature will be measured by the following shape function: 

Definition 3.2 Let (y9 : M — ]R_|. be a nonnegative even function (f{—k) = (f{k) > for 
k ^ such that ip{k) is nondecreasing for k > 0. 

As an example of a shape function one can consider (p{k) = 1, (f{k) = \k\, ip{k) = 
1 — e + e\k\ {e G [0,1]), ip{k) = \/e^+k^ {\e\ ^ 1), etc. In numerical computations 
contained in this paper, a linear combination of these functions will be often used (see 
Example 13. 7p . 

In order to answer Question 13.11 we introduce a generalized relative local length 
adopted to the shape function as follows: 

r>,t) = -^^^^^^, ..€[0,1], tG[0,T). (3.1) 

Here the bracket function (F) denotes the average of function F{u,t) along the curve F*: 
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Next we explain the role of the generalized ratio r<^. Suppose, for a moment, that r(^(u, t) = 
1 for all u at a time t. Then s{ui,t) — s{ui-i,t) ^ V^/N holds, if \k\ satisfies '^{k) ^ {^{k)) 
on the 2-th interval [ui^i,Ui]. Using this property we are in a position to construct the 
desired curvature adjusted tangential redistribution which takes account of deviations of 
the shape function ip{k) from its average value {ip{k)). Indeed, suppose that we are able 
to construct a tangential velocity a in such a way that r<^(-u,t) — )■ 1, i.e., 

9{u,t) = \nr^{u,t) 

holds for all u G [0, 1] as t — > T. For t close to T, we can conclude if \k\ is above/below 
the average in the sense that ip{k) ^ {(p{k)), then g ^ holds, respectively. Distribution 
of grid points on corresponding subarcs is dense/sparse, respectively. 

Example 3.3 Consider (p{k) = 1. Then 6' = Inr. If one constructs a in such a way that 
6{u, t) = 6{u, 0) for all u and t then relative local length r is preserved. On the other hand, 
if one takes a such that 6{u,t) — (i.e., r — )■ 1) for all m as t — > T then redistribution 
becomes asymptotically uniform. See Example 13.51 

Convergence limt^T '^ip{u,t) — )■ 1 is fulfilled provided that there exists a relaxation 
function to E Ll^^[0,T) such that 

T 

uj{t) dr = oo 

and e{u,t) = ln((e^("'°) - l)e~i'<NM°'^ + 1) for all u G [0, 1] and t G [0,T). The previous 
equation can be rewritten as an ODE: 

dte{u,t)+uj{t){l-e'^) = 0. (3.2) 



Remark 3.4 For 6 = \nr^ one can obtain another convergence such as 6{u, t) = 6{u, 0)e~ '^('^) 
then we have the equivalent ODE dt9{u, t)+uj{t)6{u, t) = 0. The convergence rate of 9{u, t) 
is the same as that of 9{u,t). 

With regard to [16l [17] one can choose the relaxation function u{t) as follows: 

u{t) = Ki — K2dt In L* = Ki + K2{k(3), 

where > and ^2 > are nonnegative constants. Here we have employed the total 
length equation: + j^t kf3ds = (cf. [15]). If we formally set Ki = K2 = 0, then 

the function 6 is constant in time t. On the other hand, if the maximal existence time is 
infinite T = oo, we can choose ki > and ^2 = 0. If T < oo and L* — )■ as t — i- T, we 
can choose K2 > 0. In both cases we have ui{t) dr = 00 and so r^{u,t) — )■ 1 as t — )■ T. 
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The equation for the tangential velocity a can be derived as follows. Differentiating 
^ = In r<^ = In (yf + In — In jj^ ipg du with respect to t, using (12. 3p . (12. 4p and the relation 

ds{fa) = a{dsk)ip' + {dsa)ip, we obtain 

^1 

y^dte = {-kp + dso)^ + {dl(5 + ei3 + adsk)^' - / dt{vg) du 

L\^) Jo 

= 9,(^a)-/ + ^(/), 

where / = ip{k)kl3 - ip'{k) {dlP + ^P), ^\k) = dkv{k). Hence holds if and only if 
the tangential velocity a satisfies the following equation: 

Equation (13. 3 p is written in the form ds{y:>{k)a) = J-", where 

^ := / - (^^(^) + ^itMm^^p' - !)• 

Since J^t — 1) ds = {L^ {ip{k)) g^^ — (p{k))gdu = 0, we obtain (J^) = 0. Thus 

the equation ds{'^{k)a) = T has at least one solution a. In order to construct a unique 
solution a, we assume the following renormalization condition for a: 

{^{k)a) = 0. (3.4) 

Example 3.5 (uniform redistribution) In the case where <f{k) = 1, the tangential 
velocity equation becomes 

d^a = k/3- {k/3) + uj{t){r^ - 1). 

If we set a; = 0, then a solution a preserves the relative local length r (see [15]). Un- 
der the assumption u ^ with suitable k,i,K2, redistribution of grid points becomes 
asymptotically uniform [161 [13 H] • 

Example 3.6 (crystalline tangential velocity) Suppose that the evolving curve F* is 
convex. If we consider the shape function (f{k) = \k\ and u{t) = 0, then, with regard to 
(13. 3p we have ds{ka) = —dlP. Taking into account renormalization constraint {(f{k)a) = 
we end up with a = —dgP/k. This is exactly the same tangential velocity as it was 
derived by the second author in the continuous limit of the crystalline curvature flow 
(see [29]). If the evolving curve F* has negative curvature on those parts, the shape 
function (p{k) = \k\ is regularized such as, for instance, (p{k) = \fl^^T^r'eW with e G 
(0,1). 
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Example 3.7 (their linear combination) In our numerical computations in SectionlHl 
we will use the following smoothed shape function ip: 



ip{k) = !-£ + ey/l - e + eP, £€(0,1). 

It is a linear combination of shape functions in Examples [33] and [321 Notice that (p{k) — t- 1 
if e ^ 0, ip{k) \k\ if e ^ 1, and ip{k) > (f{0) > for e G (0, 1). 

In Figure 13.11 we plot an example of distribution of grid points along the ellipse. 
Let x{l) = (a cos(2'7r/), 6sin(27rZ))'^ be the ellipse with axes ratio a = 3 : 6 = 1, where 
/ G [0, 1] C M/Z. We will construct a new parameterization u G [0, 1] using the reparame- 
terization function l{u) such that x(u) = x{l{u)), /(O) = 0, /(I) = 1, dul{u) > as follows: 
since dux{u) = dix{l{u))dul{u), we have g{u) = g{l{u))dul{u), where g{l) = \dix{l)\. Then 
if r<^ = 1, we obtain the ODE 



to be solved. We have used the relation k{u) = k{l{u)). Applying the Runge-Kutta ODE 
solver we calculated results depicted in Figure WA\ (a), (b) and (c). Crystalline curvature 
redistribution {x{li)} is obtained from the condition t{li) = (— sin(27ri/A^), cos(27rz/A^))'^, 
and the z-th inward normal vector of the circumscribed polygon is n(/j) for i = 1, 2, ■ ■ ■ , 
(Figure [311(d)). 



4 Optimal redistribution of points for stationary curves in the plane 

The aim of this section is to further motivate the study of curvature adjusted tangential 
redistribution. In what follows, we will address the question on what is the optimal 
redistribution of vertices {xi,X2, ■ ■ ■ ,xn} belonging to a given smooth closed curve F 
such that the discrepancy between the length/the area of F and that of a polygon spanned 
by vertices {xi, X2, ■ ■ ■ , x^} is minimal. Clearly, in the case F being a circle, the optimal 
redistribution is represented by a regular A^-polygon. However, it should be obvious 
that, e.g. for an oval, the optimal redistribution has to take curvature of the curve 
into account. Interestingly enough, we will prove that the length/the area discrepancy 
minimizing redistribution X = {xi,X2, ■ ■ ■ ,xn}, for — )■ oo, is closely related to the 
curvature adjusted redistribution discussed in previous section with the shape functions 

tf{k) = m"^^^ for the length discrepancy, 
ip{k) = \k\^^^ for the area discrepancy. 

Suppose that a given smooth closed curve F is parameterized by a smooth function x : 
[0,1] I— )■ M}. Denote by X = {xi,X2,--- ,X]\f} the set of grid points of F such that 
Xi = x[ui) where Ui = ih and h = 1/N. 
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The length discrepancy. Let L = J-pds be the length of F and 

be the length of a polygon with vertices {xi, X2, ■ ■ ■ , x^}. Here we have identified x^^i = 
Xi since F is assumed to be a closed curve. Clearly, C{X) < L. Our goal is to find condi- 
tions under which the parameterization x[-) yields the minimizer X = {xi, X2, ■ ■ ■ , x^} 
of the problem 

min(L-/:(X)). (4.1) 

A Ci 

Since 

i-2 N-1 

we obtain the following expression for the derivative of C with respect to Xi in the direction 
y e M^: 

a^AXjy = -rii.y, rii 



Recall that the above minimization problem (14. ip is subject to the constraint a?j G F for 
each i = 1, 2, ■ ■ ■ ,N. Therefore, the first order necessary condition for the constrained 
minimizer X of (14. ip reads as follows: 

C'^^{X)y = foranyt/~t„ z = 1, 2, ■ ■ ■ , iV, 

where the symbol y ^ ti means that the vector y is collinear with the unit tangent vector 
ti to F at the point Xi G F. Hence, the necessary condition for a set X = {a^i, ■ ■ ■ ^x^} G 
]^2xAr q£ points belonging to F to be a minimizer of the functional X L — C{X) can be 
rewritten as: 

Ui^ti, z = 1,2,--- ,Ar. (4.2) 

A graphical description of the necessary condition is depicted in Figure WA\ (a). 

Next we will express the necessary condition (14. 2 p in terms of the parameterization 
£c(-) of the curve F. The Taylor expansion of a:;(-) at Xi = x{ih) yields 

fi'i fi-i 

Xi±i -Xi = ±dux{ui)h + dlx{ui)— ± dlx{ui)— + d^x{ui)— + 0{h^), 

Z b 24: 

as — 7- 0^. Derivatives d^x at Ui = ih can be decomposed as follows: 
dix{ui) = bjUi + Qjti, j = 1, 2, ■ ■ ■ , 
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where and tj are, respectively, the unit inward normal and tangent vector to the curve 
r at a point Xi G F. Then ai = 61 = where g = Qi = \dux{ui)\ is the local length 
at the point Xi G F. Furthermore, using the Frenet formulae g~^dut = dgt = kn and 
g^^duTi = dsTi = —kt, we obtain the recurrent relations 

ai = g, bi = 0, aj_|_i = duttj — gkbj, bj+i = dubj + gkaj. (4.3) 

for j > 1 where k = ki is the curvature of F at ajj G F. 

Let us define and "^{h) as the following auxiliary functions 

00 °° I, 

i=i i=i 
Using these functions, the forward and backward difference at Xi can be expressed as 



Then the necessary condition (14 .21) can be rewritten as: 
F{h) = 0, 

where h = 1/N. The function F is defined by 
F(h) = '^(^) _ '^(-^) 

^$(/l)2 + ^(/l)2 ^$(-/l)2 + ^(-/i)2' 

Due to the symmetry F{-h) = -F{h) we have F(0) = F"(0) = 0. Moreover, F'(0) = 
because of the fact that 61 = and so ^E'(O) = 0. The leading order term of Taylor's 
expansion of F{h) is therefore given as: 

F{h) = ^F'"(0)/i=^ + 0{h^) 
o 

as /i — )■ 0+. Calculation of F"'(0) can be simplified if one uses the substitution ^(h) : = 
\&(/i)/<l>(/i). Notice that ^{h) > for small < 1 because $(0) = ai = > 0. Then 
F"'{0) = -6^'(0)^"(0). Using we finally obtain 

F"\0) = 3^ [- - 

^ ' all 2 3 ^2 

Since F{h) = we finally deduce the condition 

02 _ 1 ai&3 _ Q 
2 3 62 ~ ' 

which has to be satisfied by the minimizer X in the limit h 0^, i.e. for — )• 00. Now 
calculating the corresponding terms ai, 02, &2 and 63 from the recurrent relations (14.31) we 
obtain 

ai = g, 0,2 = dug, &2 = %^ &3 = du{kg'^) + gkdug. 



11 



Hence 

Id^g 1 du{kg^) + gkdug 



on r = {x{u); u E [0, 1]}, 



2 g 3 kg"^ 
which is clearly equivalent to the statement: \k\g^^^ = constant, i.e. 

I A; 1 2/3^ = constant (4.4) 

on the curve T = {x{u); u G [0, 1]}. It corresponds to the curvature adjusted tangential 
velocity with '^{k) = and the extended relative local length (see fl3.ll) ) 

r,(n) = ^f^^ = 1, for each [0,1]. 
L (<^(fc(-))) 

In other words, such a curvature adjusted tangential velocity yields (in the limit of number 
of grid points tending to infinity) the best possible approximation of the evolving family 
of planar curves as far as the minimization of the error between the length of a curve and 
the length of its polygonal approximation is concerned. 

Figure [XT] (e) displays an example of the length discrepancy minimizing redistribution 
for the case when 

^(k) = |A;|2/3. 

The values of length defect Al := 1 — C{X)/L for several X's corresponding the cases in 
Figure 13.11 and the length optimal case are calculated as in Table 14.11 

The area discrepancy. Similarly as in the case of the length functional we can 
ask the question what is an asymptotically optimal redistribution of points on a given 
curve r such that the discrepancy in areas enclosed by the curve F and its polygonal 
approximation is minimal. 

The area A enclosed by the curve F is given by A = | det{x, dgx) ds. The area 
enclosed by the closed polygonal curve with vertices X = {xi, X2, ■ ■ ■ , a^Ar} is given by 

= 2 5Z d6t(a;j, Xj+i - Xj), 

3=0 

where Xq = x^. Clearly, A{X) = \J2f=o ^^'^i^j^^j+i)- The first order necessary condi- 
tion for a set X = {xi, X2, ■ ■ ■ , a; at} of vertices to be a minimizer of the area discrepancy 
functional 



mm(A - A(X)) 

xcr 



where the vertices are constrained to belong to the curve F reads as follows: 
A'^^iX)y = ioianjyr^t,, ^ = 1, 2, ■ ■ ■ , iV. 
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We have 

,1 1 _ _ 

A'^^{X)y = - {det{xi^i,y) + det{y, Xi+i)) = --det{ti,y), U = Xi+i - Xi^i 

for any direction y G M^. Hence, the necessary condition for a set X = {xi, ■ ■ ■ ,xj\f} G 
■[^2xAf q£ points belonging to F to be a minimizer of the functional X {A — ^(X))^ can 
be rewritten as: 

ti II U ^ ti±ni, t = l,2,--- ,N. (4.5) 

A graphical description of the necessary condition is depicted in Figure W7[\ (b). 

Now, calculating the difference ti by means of the Taylor series and using the facts: 
det(ti,nj) = 1 and det(ti,ti) = we obtain that condition (14.51) can be restated as: 

h,h^ + oih'') = o 

as /i — > 0^. Therefore, in the limit h = l/X — )■ 0, we conclude that 63 = 0. Since 
^3 = du{kg^) + gkduQ = kg'^du{\'^{\k\g'^) H-ln^f), we obtain g'^\k\ = constant on F, i.e. 
5f|/i;|^/^ = constant. It means that the area discrepancy curvature adjusted tangential 
velocity has the shape function 

ip{k) = \k\^/\ 

Figure Em (f) displays an example of the area discrepancy minimizing redistribution 
for the case when (fi^k) = \k\^^^. The values of area defect := 1 — A{X)IA for several 
X's corresponding the cases in Figure [3T] and the area optimal case are calculated as in 
Table O 

5 Numerical approximation scheme 

The purpose of this section is to construct a numerical approximation scheme for solving 
the governing equation (12. 2p for the position vector and the tangential velocity equation 
(13. 3p satisfying renormalization constraint (13. 4p . 

Scheme. For a given initial X-sided polygonal curve 

^° = UiLi find a family of 
iV-sided polygonal curves {'P-'}j=i,2,--, = [JiLiSj, where S{ = [xi_-^,xi] is the i-th 
edge with x^ = x^'jy for j = 0, 1, 2, ■ ■ ■ . The initial polygon is an approximation of 
F° satisfying {x^ijfL^ C P° n F°, and is an approximation of F* at the time t = tj, 
where tj = jr is the j-th discrete time (j = 0, 1, 2, ■ ■ ■ ) if we use a fixed time increment 
r > 0, or tj = Xlto (j = 2, ■ ■ ■ ; to = 0) if we use adaptive time increments ti > 
0, / = 0, ■ ■ ■ , j — 1. The updated curve {V^^^} is determined from the data {V^} at the 
previous time step by using discretization in space and time as follows. 

Discretization in space. Let V = IJili '^j t)e an iV-sided polygonal curve, where 
Si = [xi_i,Xi] is the i-th edge and Xi is the z-th vertex (z = 1,2, ■■■ ,A^; Xq = x^). 
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We denote the length of Si by = \xi — Xi_i\. The i-th unit tangent vector can 
be defined as U = {xi — Xi^i)/ri. Then the i-th unit tangent angle z/j is obtained from 
ti = (cos z/j, sin z/j)"'" in the following way: firstly, from ti = (tii,t2i), we obtain i>i = 
2tt — arccos(tii) if < 0; ui = arccos(tii) if ti2 > 0. Secondly, for z = 1, 2, ■ ■ ■ , we 
successively compute z/j+i from uf. 

{Vi + arcsin(D) if / > 0, 
Vi + arccos(/) if D > 0, where D = det(ti, t^+i), / = tj-tj+i. 
Ui — arccos(/) otherwise. 

Finally, we obtain uq = ui — {i^n+i — ^n) and z^Ar+2 = t^n+i + {^2 — ^i)- 

In order to derive a discrete numerical scheme, we follow the flowing finite volume 
method adopted for curve evolutionary problems as it was proposed by Mikula et al. 
in [T71 [IS]. Let us introduce the "dual" volume S* = [x*,Xi] U [a^j, a;*_,_^] of Si, where 
X* = {xi_i + a3j)/2 (i = 1, 2, ■ ■ ■ , A^; x*j^j^^ = x\). We define the z-th unit tangent angle 
of S* by V* = {ui + z/j+i)/2. The z-th curvature ki has the constant value on Si, which is 
obtained from integration oi k = dgV over Si with respect to s: 



I kds = ki ds = kiTi, I kds= I dgV ds = [i'Yxi_^ 

J Si ^ Si ^ Si ^ Si 



i-l- 



Hereafter, F ds means J^'^' ^ F ds for arc-length Si satisfying Xi = x{si, ■). Thus we have 
ki = {dsi'*)i, where {dsf)i = (Fj — Fj_i)/rj. The i-th curvature k* at Xi can be defined as 
k* = {ki+i + ki)/2 which has the constant value on S*. 

Next we discretize equation (I3.3P for the tangential velocity a: 

where / = (9^/3 + k'^(3)ip'{k) — k(3ip{k). Integrating the above equation over Si yields 

ds{ipa)ds = [^{k)a]l'^_^ 

= V{k)ds-j fds+(^L{^) j ^ds- [ ^{k)ds^u. 



Then 



ipi = ^{k*)ai - <^(A;*_i)a 



i-l 



^■^'^ ip{ki)ri - fiVi + ( C{ip)^ - (p{ki)ri ] co, 



where Pi = /3{x*, ki,ti) is constant on Si, 
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r* = {vi + ri^i)/2 is the length of S*, and 

N 



N 
i=l 



is the total length of V. The averages (/) and (ip) are approximated as: 

1 ^ 1 ^ 

i=l i=l 

To determine we have to take account of the renormalization constraint 

(fa) = 0. It can be discretized as: 

1 ^ 
^ i=i 

Notice that £ = J2iLi = Xlili ^i- We define a partial sum of {ipi} by 

i 

^, = 5^^'; (^ = 2,3,--- ,iV), *i = 0. 
With this notation we obtain 

Summing the above terms yields 

N N N 

i=l 

Hence we obtain 



i=l i=l 1=1 



N 

i 



1 ^ 1 



i=2 



Discretization in time. The semidiscretized (in space) evolution equation fl2.2p has 
the form: 

dtXi = w*{ds*t)i + ai{ds*x*)i + F*n*, w* = w{xi, z/*, k*) F* = F{xi, v*) 

for z = 1, 2, ■ ■ ■ , A^, which follows from integration of dtX = wdgt + adgX + Fn (see (12.21) ) 
over the volume S*. The right-hand side can be discretized as follows: 

9.^. - ^ (^^ - + - -.0 + F-m. 
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In our approach, we use the semi-imphcit numerical scheme for discretization in time: 



— — a, do = — (a_ + a+) 



where is the i-th data of an iV-sided polygonal curve ~ ^{tj) for ^ = 1, 2, ■ ■ ■ , A^. In 
other words, the following tridiagonal matrix has to be solved in order to obtain solution 
in the new j + I time level: 

-a^Txi^l + (1 - aor)xf ^ - a+rx^^^ = xl + F^^ul^r = 1, 2, ■ ■ ■ , iV) 

subject to periodic boundary conditions xi^^ = x-^j^^, a^i^^ = x{^^. We note that the 
tridiagonal matrix is strictly diagonally dominant provided that r is small enough. The 
time step r can be also chosen adaptively using the following expression 

where A > 0, Fmin = mini<i<Ar Fj, F^ax = maxi<i<Ar F^ and |F|max = maxi<i<Ar |Fi|. With 
this choice of tj the solution {a;^^^}^^ exists and it is unique. 



6 Computational results 

This section is devoted to presentation of various numerical experiments. In all figures, for 
the prescribed r > 0, we plot every /ir discrete time step using discrete points representing 
the evolving curve. In every 3yur time step, we plot a polygonal curve connecting those 
points, where n = [[T/r] /lOO] {[x] is the integer part of x), and T is approximation of the 
final time which computed as follows. Let and £* be the enclosed area and the total 
length of at time t = tj. In Figures [6.71 16.81 and 16. 9^ T is the smallest discrete time t 
for which both conditions \A'^/A^^'^ — 1\ < 6 and \0' j C}^'^ — \ \ < b are satisfied. In all 
other figures, T is the smallest discrete time t such that j < 6. In figure captions, 
we provide the number of grid points N, e G (0, 1) for (p{k) = 1 — e + eVl ~ s + sk?, ki 
and K2 for the relaxation function uj{t). In all examples, the initial discretization is given 
by the uniform A^-division of the u-range [0, 1]. 

Initial test curves. As an initial test examples we use a circle, an ellipse and the 
following initial curves x{u, 0) = {xi{u) , X2{u)y- depicted in Figure [6TT] and parameterized 
by 

Xiiu) = cosz, X2{u) = 0.7sin2; + sina;i + Xg, u E [0, 1], 
where x^ = sin(3z) sinz and z = 27tu (left), and 

xi{u) = 1.5 cosz, X2{u) = 1.5(0.6 sin z + 0.5x1 + 0.4 sin 0:4 + 0.1 sin 0:5) , u E [0, 1], 
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where 0:3 = sin(32;) sin 2;, X4 = 2x1, ^5 — 3e~^i and z = 27ru (right) (cf. [T3]). 

Classical curvature flows. According to the classical convexification theory and the 
asymptotic behavior derived by M. A. Grayson, M. Gage and R. S. Hamilton, in the case 
where the geometric equation is given hj P = k, any embedded curve becomes convex in 
finite time [9], and any convex curve shrinks to a single point. Its asymptotic shape is a 
circle [S] . See Figure 16.21 

Afline curvature flows. Convexification also holds true in the case where the evo- 
lution law is the so-called affine scale space normal velocity (3 = k^^^ (cf. |2l])- In this 
case, convexification was proved and the asymptotic behavior was derived by G. Sapiro 
and A. Tannenbaum. They showed that any embedded curve shrinks to a single point 
with an ellipse as the asymptotic shape [21] . See Figure 16.31 

Experimental order of convergence (EOC). In the case of the normal velocity 
P = k^^^, any ellipse shrinks to a point homothetically. It means that an ellipse is a 
self-similar solution to ( II. ip (see Figure El]). By using the explicit self-similar solution, 
the so-called experimental order of convergence (EOC) can be obtained in the following 
way. 

When the position vector is described by x{u,t) = ri[t)z{u) with ri{0) = 1, and when 
the initial curve is an ellipse x{u,0) = z{u) = {acos{27cu),bsm{2TTu))^ with a,b> 0, for 
the curvature we obtain k{u,t) = abr]{t)^^((u)^^^'^, ((u) = sin^(27r'u) + 6^ cos^(27r'u). 
In the case when (3 = k^^^, dtX ■ n = {aby^^r]{t)^^^^(^(u)^^^'^ holds. On the other hand, 
we have dtX ■ n = dtri{t) z ■ n = —abdti]{t) C{u)^^^'^. Hence we obtain the rate r]{t) = 
(1 — |(a6)~^/^t)^/^ and the extinction time T = |(a6)^/'^, which is determined by ri(T) = 0. 

Numerical errors at t = tj with the number of points can be defined as 



max 

l<i<N 



eiiUN) = < 



[X- 



J \2 



+ 



[X. 



J \2 
2i) 



- y 



- 1 



\<i<N 



(ar/(t,))2 {b^{t,)f 



p\ i/p 



if p = 00, 



if 1 < p < 00, 



where xj = {x-' ^ 
error such as 



^ii^^2i)" is the z-th grid point. Therefore we can define the L«((0, Uf), ^^^(0, 1)) 



max err 



if g = 00, 



1/9 



— J2 i'^Ky] ' ifl<g<oo. 



i<i<Af 



where Im < T. If Ep^g{N) ^ N~^' holds, then Ep^g{N/2) ^ 2^'N-f also holds. From 
these approximations, we may expect /i ~ logg [Ep^g(N / 2) / Ep^g^N)). We therefore use 
the right hand side as an indicator of the numerical convergence, i.e., the EOC from the 
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Li{{0,tM),LP{0, 1)) error is defined as 



EOCp,,(iV) 



E 

log 



{N/2) 



Tables EB E21 E31 ESI indicate the L«((0, tM), ^^^(0, 1)) errors ^p,g(A^) and EOC 
EOCp,g(A^) for N = 2^ 2^ 2^, 2\ 2^ (p, g) G {1, 2, oo}, and e = 0, 0.1, 0.5, 0.9. Parameters 
are r = 0.1A^~^ and ni = k,2 = 100, and sample data are used at the time tj = 1.5j/M 
with j = 0,l,---,M = 200. Ellipses with the axes ratio 3 : 1 are the same as in Figure |6]H 
In this case, the extinction time T = |(a6)^/^ = 1.560 ■ ■ ■ > = 1-5. Since the governing 
system is of the parabolic type, we have used the time step r ^ N'"^. From Tables 16. 1[ 
16. 2[ 16. 3[ 16. 4[ we can observe EOC ~ 2 independently from the values e's. It means that 
these four tangential redistribution methods seem to be almost equally effective in view 
of the EOC. 

The length and the area discrepancy. An ellipse is a shrinking self-similar solution 
to (11. ip in the case of the normal velocity (3 = w{h')k, as well as the affine curvature flows 
(3 = k}^^. Here the weight is 



w(z/) 



2T{a? sin^ z/ + 6^ cos^ p) 



when the axes ratio of ellipse is a : 6 and the extinction time is T > 0. Similarly as in 
the previous subsection EOC, one can obtain the position vector of a shrinking ellipse by 
x{u,t) = rj{t)z{u) with the scaling function rj{t) = a/1 — t/T (77(0) = 1, ri{T) = 0) and 
the initial ellipse x{u, 0) = z{u) = {a cos{27iu) , 6sin(27rM))^. By using this exact solution, 
we will calculate the numerical test of the length and the area discrepancy in the following 



way. 



Numerical discrepancy of the length and area defects at the time t = tj can be defined 



as 



1 



A-' 



1 



Ah 



respectively. Here L* = r]{t)L^ and A* = ri{t)'^A^ are the length and area of F*, and 
and A* are the length and area of at t = tj, respectively. Therefore we can define the 
L'^{0,tM) numerical discrepancy of the length and follows 



max Ai; 




if q 



00, 



1/9 



0<j<M 



if 1 < g < cxD, 



where * = L,A and tM T. 

Table 16.51 indicates A and A^ ^ for q 



1, 2, 00. Parameters are: = 100, r 



O.IA^ ^ and ki = k,2 = 100, the extinction time is T = 1, sample data are used at the 
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time tj = Tj/M with j = 0, 1, ■ ■ ■ , M = 200, and the axes ratio of the initial elhpse r° 
are a = 3 : 6 = 1. As for the shape functions we chose: v9(/c) = 1 — e + e\J\ — e + ek"^ 
with e = 0,0.5,0.9, 1, ^(k) = \k\'^^^ and ^(k) = \k\^^^. In Table [675| we can observe that 
the shape function (p{k) = \k\ attains the minimum value in each A*^g for * = L,A and 
q = 1,2, oo. 

Weighted curvature flows. Asymptotic behavior of solutions to the weighted cur- 
vature flow P = w{v)k is related to self-similar shrinking solutions, which need not be 
unique. For details we refer to the paper by Yagisita [28] and references therein. In 
Figure 16751 we show evolving family of plane curves with /3 = w{i>)k. 

Forced curvature flows. The theory of Grayson, Gage and Hamilton was generalized 
to the case of an anisotropic curvature driven flow by K. S. Chou and X. P. Zhu in [3]. 
They considered the case where the evolution law is given by /3 = w{v)k + F{v) with 
w(z/) = a{u) + dlaiv), where a is a given anisotropy function, cr(z/ + tt) = cr(i^), and 
F{v + tt) = —F{v). They proved that any embedded curve becomes convex in flnite 
time [4], and any convex curve shrinks to a single point with the asymptotic shape being 
a self-similar solution to (3 = w{i')k, i.e., the Wulff shape of a [3J. In Figure [676] we 
illustrate the convexification theory and the asymptotic behavior of evolving plane curves 
due to K. S. Chou and X. P. Zhu. 

Loss of convexity phenomenon. In the case when the normal velocity is given 
hy (5 = k + F{x, v), an initially convex curve may loose its convexity in flnite time for 
a special choice of the external force F. In Figure 16.71 we present such examples with 
F presented by Nakamura et al. in [21]. It is worth noting that the usual crystalline 
curvature flow is unable to capture this convexity-breaking phenomena. 

Image segmentation by using a gradient flow. Let 7(3;) > be an inhomogeneous 
energy density along the curve F*. If 7 is differentiable, then the gradient flow of the 
following energy 

E{T') = [ -f{x)ds 

is realized by the geometric equation (11.11) of the form /3 = '~f{x)k — V'y{x) ■ n. Such a 
gradient flow can be successfully utilized in various image segmentation problems. For 
example, let an image intensity function be denoted by / : D ^2 — )■ [0, 1]. We assume 
that I{x) is a piecewise constant function on each pixel. Here J = (/ = 1) corresponds to 
the black (white) color and / G (0, 1) corresponds to a scale of gray colors. For simplicity, 
we assume that our target flgures are drawn in white color with the black background. 
Then edges of the image correspond to regions where the gradient |V/(a3)| is sufficiently 
large. Let us introduce an auxihary density function 7(03) = f{\VI{x)\) where / is a 
smooth edge detector function such as /(s) = 1/(1 + s^) or /(s) = e~'^. Hence the 
solution curve F* of /3 = '~f{x)k — V7(a;) ■ n has the tendency to minimize the energy 
i?(F*). In other words, it moves toward the edge in the image on which |V/(a3)| is large. 
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This is a fundamental idea of image segmentation, and it has developed to a sophisticated 
numerical scheme [TTJ [18]. An example of such a curve evolution converging to an edge 
in the given images is depicted in Figure 16.81 

Image segmentation for sharp images. If contrast of the target image is sufficiently 
high, a simpler scheme described in |2j can be used. We consider the geometric flow 
P = k + F and define the forcing term F{x) as follows: 

-^{^•^) -^max (-^max -^min) -^(''^) ("^ ^ 

where F^ax > corresponds to pure black (background) and F^i^ < corresponds to 
pure white (the object to be segmented). Maximal and minimal values determine the 
final shape because in general 1/F is equivalent to the minimal radius the curve can 
attain. Choosing an appropriate F, we can make the final shape be rounded, or we can 
prevent the curve from passing through the outline. A segmentation of a given sharp 
image by means of a plane curve evolution is shown in Figure 16.91 

7 Conclusions 

In this paper, we proposed and analyzed a new class of tangential velocities by which we 
can control tangential motion and grid point redistribution of plane curves that evolve with 
the normal velocity depending on a general function of the curvature, tangential angle and 
the position vector. The curvature adjusted tangential velocity may not only distribute 
grid points uniformly along the curve but also produce a desirable concentration and/or 
dispersion depending on the curvature. We also demonstrated that curvature adjusted 
tangential redistribution yields the best possible approximation of the evolving family 
of planar curves as far as the minimization of the error between the length and area 
of a stationary curve and that of its polygonal approximation is concerned. Numerical 
experiments based on semi-implicit numerical fiowing finite volume method confirmed 
capability of our new method. 

Acknowledgments. The authors would like to express our gratitude to the anonymous 
referees for their comments and suggestions. 
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(e) (f) 

Figure 3.1: Various redistributions of = 12 grid points along the ellipse: (a) uniform re- 
distribution (e = 0), (b) curvature adjusted redistribution {e = 0.9), (c) curvature adjusted 
redistribution (e = 1), (d) crystalline redistribution (corresponding to e = 1), (e) the length 
discrepancy minimizing curvature adjusted redistribution with (p{k) = \k\'^^^ , and (f) the area 
discrepancy minimizing curvature adjusted redistribution with (p{k) = \k\^/^. 



X 


(a) uniform 
e = 


(b) curv. adj. 
e = 0.9 


(c) crystalline 
e = 1 


(e) the length 
optimal 


(f) the area 
optimal 


Al 


0.01835 


0.00789 


0.00966 


0.00733 


0.01085 




0.05834 


0.05400 


0.11998 


0.06222 


0.04507 



Table 4.1: The values of length defect Al and area defect for the cases in Figure [3T] (a), 
(b), (c), (e) the length optimal case ip{k) = |A;|2/3, and (f) the area optimal case (p{k) = \k\-^''^. 
Underlined values are the minimum in each defect. 
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Figure 4.1: A graphical illustrations of the necessary conditions at a point Xi for maximizing 
of the polygonal length discrepancy C{X) (a) and the polygonal enclosed area A{X) (b) of a 
discrete approximation X = {xi, - ■ ■ ,xj\f} of a smooth curve F, respectively. 



Figure 6.1: Initial curves with grid points corresponding to a uniform division of u G [0, 1] with 
N = 100 points. 




(a) (b) (c) (d) 

Figure 6.2: Numerical examples with the normal velocity (3 = k. (a) and (c) are evolving curves, 
respectively, (b) and (d) are the corresponding final magnified curves. Numerical parameters: 
N = 100, e = 0.1, Ki = K2 = 100, r = O.ITV'^, ,5 = lOr, and f = 0.001. 
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(a) (b) (c) (d) 

Figure 6.3: Numerical examples with affine the space scale normal velocity /? = k}/^: (a) and 
(c) are evolving curves, and (b) and (d) are the corresponding final magnified curves. In both 
cases, parameters were chosen: N = 100, £ = 0.1, ki = K2 = 100, r = O.IA*""^, 5 = lOr, and 
f = 0.001. 




Figure 6.4: From (far left) to (far right), e = 0,0.1,0.5,0.9. Upper figures are evolving ellipses 
starting from the same ellipse with the axes ratio 3:1. Lower figures are the final magnified 
curves. Numerical parameters: N = 128, k,i = K2 = 100, r = 0.1A^~^, 5 = lOr, and r = 0.001. 
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N 


P 




EOCp,i{N) 




E0Cp,2(iV) 




EOCp,oo(iV) 


16 


1 


0.5380220 




0.5509331 




0.6336300 






2 


1.0534598 




1.0704854 




1.2142337 






oo 


5.9554974 




6.0404891 




6.8679124 




32 


1 


0.1792467 


1.586 


0.1825282 


1.594 


0.2122894 


1.578 




2 


0.3427337 


1.620 


0.3437584 


1.639 


0.3692181 


1.718 




oo 


1.8652997 


1.675 


1.8668395 


1.694 


1.9453397 


1.820 


64 


1 


0.0537848 


1.737 


0.0552922 


1.723 


0.0691816 


1.618 




2 


0.1018371 


1.751 


0.1024454 


1.747 


0.1161671 


1.668 




oo 


0.5383176 


1.793 


0.5391202 


1.792 


0.5812710 


1.743 


128 


1 


0.0148900 


1.853 


0.0154231 


1.842 


0.0202166 


1.775 




2 


0.0281297 


1.856 


0.0283963 


1.851 


0.0334721 


1.795 




oo 


0.1466961 


1.876 


0.1470527 


1.874 


0.1622977 


1.841 


256 


1 


0.0038891 


1.937 


0.0040411 


1.932 


0.0053923 


1.907 




2 


0.0073462 


1.937 


0.0074271 


1.935 


0.0088928 


1.912 




oc 


0.0381101 


1.913 


0.0382-141 


1.913 


0.01255C0 


1.931 



Table 6.1: e = 



16 


1 


0.5290974 




0.5419132 




0.6255006 






2 


1.0325342 




1.0494491 




1.1963774 






oo 


5.8214106 




5.9054873 




6.7533333 




32 


1 


0.1743048 


1.602 


0.1774370 


1.611 


0.2056751 


1.605 




2 


0.3314933 


1.639 


0.3324330 


1.658 


0.3553780 


1.751 




oc 


1.7951105 


1.G97 


1.7901082 


1.717 


1.8597097 


1.801 


64 


1 


0.0518231 


1.750 


0.0532112 


1.738 


0.0656620 


1.647 




2 


0.0974615 


1.766 


0.0979813 


1.762 


0.1092110 


1.702 




oo 


0.5117050 


1.811 


0.5122756 


1.810 


0.5425148 


1.777 


128 


1 


0.0142197 


1.866 


0.0146984 


1.856 


0.0188710 


1.799 




2 


0.0266349 


1.872 


0.0268569 


1.867 


0.0309018 


1.821 




oo 


0.1377070 


1.894 


0.1379516 


1.893 


0.1485155 


1.869 


256 


1 


0.0036908 


1.946 


0.0038251 


1.942 


0.0049766 


1.923 




2 


0.0069038 


1.948 


0.0069696 


1.946 


0.0081078 


1.930 




oo 


0.0354977 


1.956 


0.0355656 


1.956 


0.0384550 


1.949 



Table 6.2: e = 0.1 
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N 


P 




EOCp,i{N) 




E0Cp,2(iV) 




EOCp,oo(iV) 


16 


1 


0.4042607 




0.4129555 




0.4910215 






2 


0.7821867 




0.7914808 




0.9293693 






oo 


4.3834563 




4.4198180 




5.1595141 




32 


1 


0.1310031 


1.626 


0.1322110 


1.643 


0.1454189 


1.756 




2 


0.2483055 


1.655 


0.2485107 


1.671 


0.2563577 


1.858 




oo 


1.3333636 


1.717 


1.3334130 


1.729 


1.3505401 


1.934 


64 


1 


0.0362174 


1.855 


0.0366554 


1.851 


0.0410071 


1.826 




2 


0.0669561 


1.891 


0.0670737 


1.889 


0.0698023 


1.877 




oo 


0.3448818 


1.951 


0.3449364 


1.951 


0.3496285 


1.950 


128 


1 


0.0093125 


1.959 


0.0094400 


1.957 


0.0106475 


1.945 




2 


0.0170431 


1.974 


0.0170821 


1.973 


0.0178686 


1.966 




oo 


0.0863334 


1.998 


0.0863494 


1.998 


0.0876178 


1.997 


256 


1 


0.0023454 


1.989 


0.0023787 


1.989 


0.0026894 


1.985 




2 


0.0042800 


1.994 


0.0042904 


1.993 


0.0044945 


1.991 




oc 


0.0215719 


2.001 


0.0215790 


2.001 


0.0218997 


2.000 



Table 6.3: e = 0.5 



16 


1 


0.4408205 




0.4417817 




0.4758601 






2 


0.8949908 




0.8959769 




0.9555951 






oo 


5.1442600 




5.1492367 




5.4758260 




32 


1 


0.1161771 


1.924 


0.1164446 


1.924 


0.1245329 


1.934 




2 


0.2229378 


2.005 


0.2230329 


2.006 


0.2323277 


2.040 




oc 


1.20G0030 


2.093 


1.20(32137 


2.091 


1.2125380 


2.110 


64 


1 


0.0290339 


2.001 


0.0290913 


2.001 


0.0307145 


2.020 




2 


0.0543111 


2.037 


0.0543233 


2.038 


0.0558738 


2.056 




oo 


0.2827757 


2.093 


0.2827892 


2.093 


0.2876229 


2.111 


128 


1 


0.0072402 


2.004 


0.0072539 


2.004 


0.0076294 


2.009 




2 


0.0134384 


2.015 


0.0134411 


2.015 


0.0137813 


2.019 




oo 


0.0690806 


2.033 


0.0690829 


2.033 


0.0700398 


2.038 


256 


1 


0.0018087 


2.001 


0.0018121 


2.001 


0.0019041 


2.002 




2 


0.0033502 


2.004 


0.0033509 


2.004 


0.0034332 


2.005 




oo 


0.0171611 


2.009 


0.0171616 


2.009 


0.0173862 


2.010 



Table 6.4: e = 0.9 



27 



^{k) 


Al,i 


Al,2 


Al,oo 


Aa,i 


Aa,2 


Aa,oo 


ipo{k) = 1 


0.014487 


0.044087 


0.461457 


0.026691 


0.099960 


1.119496 




0.013191 


0.039820 


0.417640 


0.024203 


0.089520 


1.000038 




0.005216 


0.016007 


0.172409 


0.009032 


0.033209 


0.370852 




0.001666 


0.004468 


0.047258 


0.002119 


0.008126 


0.092089 


Vi{k) = \k\ 


0.001051 


0.001213 


0.006871 


0.001291 


0.002045 


0.018642 


|/t|2/3 


0.001124 


0.002903 


0.030946 


0.001427 


0.005267 


0.060430 




0.003680 


0.011821 


0.128524 


0.006159 


0.023999 


0.270145 



Table 6.5: Numerical discrepancy of the length ^L,q and the area ^A,q for g = 1, 2, oo in various 
choice of if{k), where (pe{k) = 1 — e + eVl — e + efe^- 




(a) (b) (c) (d) 




(e) (f) (g) (h) 



Figure 6.5: (a), (b), (c), (d) are in the case where the weight is w{v) = 1 — (7/9) cos(3z^), and 
(e), (f), (g), (h) are those where w{v) = 1 — 0.8cos(4(i/ — 7r/4)). Both weights are used in |15j. 
(a) and (c) (resp. (e) and (g)) are evolving curves, respectively, (b) and (d) (resp. (f) and (h)) 
are the corresponding final magnified curves. In all cases, parameters are = 100, e = 0.1, 
Ki = K2 = 100, T = O.liV-2, 5 = lOr, and 9 = 0.001. 
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(a) (b) (c) (d) (e) 

Figure 6.6: Numerical examples for the normal velocity /3 = w{v)k + F{u)^ w{v) = 1 — 
0.7cos(6z^), F{v) = sin(z^). (a) and (c) are evolving curves, respectively, (b) and (d) are the 
corresponding final magnified curves. In both cases, parameters were chosen as: N = 100, 
e = 0.1, Ki = K2 = 100, T = 0.1Af"2^ 5 = lOr, and r = 0.001. The function a{u) = 1 + 
(0.7/35) cos(6z^) is the unique solution oi w = a + d^a, and (e) is the locus of the boundary of 
the Wulff shape of a, which is given by n) + {dua)t. 




Figure 6.7: (a) and (b) correspond to the external force F = —2pq sm{q[4:x1 +X2))(— 4j;i sinz^ + 
X2Cosz^) with p = 1.25, q = 3.0. (c) and (d) correspond to the force F = 2pqTT cos{q7:\x\'^)x • n, 
p = 1.956, q = 1.15. (a) (resp. (c)) indicates evolving curves starting from the unit circle (resp. 
ellipse with axes ratio 1 : 2), and (b) (resp. (d)) is the final curve at T. Parameters are N = 120, 
£ = 0.5, Ki = 100, K2 = 0, T = f = 0.001, and 5 = 10~^. 
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-1.5 -1 -0.5 0.5 1 1.5 -1.5 -1 -0.5 0.5 1 1.5 -1.5 -1 -0.5 0.6 1 1.5 



(a) (b) (c) 




Figure 6.8: (a) indicates the original bitmap image faded shadow "Venus" in J7 = [—1.5, 1.5]^ 
with the 600px resolution, (b) indicates image intensity function I{x) in 3D, and (c) is the gray 
scaled auxiliary function 7(2;). (d) indicates evolving curves starting from circle with radius 1.5, 
and (e) is the final curve at T. Parameters were chosen as: N = 100, e = 0.1, ki = 100, K2 = 0, 
r = 0.001, and 5 = 10^^. Adaptive time step sizing ()5.ip with A = 1 has been used. 
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Figure 6.9: (a) indicates the original bitmap image "ya" (it is the Chinese character for arrow) 
in Q, = [—1.5,1.5]^ with the 600px resolution, (b) indicates evolving curves starting from circle 
with radius 2, and (c) is the final curve at T. Parameters are = 200, e = 0.1, ki = 100, 
K2 = 0, T = 0.0005, 6 = 10"^/2, Fmax = 30 and Fmin = -30. Adaptive time step sizing (151]) 
with A = 1 has been used. 
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